Metodos Analíticos 2023

  • Vianey Galindo Añel
  • Alberto Fuentes Chavarría

Predicción inmediata para la venta de casas usando la paquetería BSTS de R y búsquedas de Google.

Contexto

La predicción inmediata (nowcasting) es una écnica poderosa que permite estimar y predecir variables económicas y sociales en tiempo real, proporcionando información valiosa para la toma de decisiones rápidas y fundamentadas. Utilizar nowcasting permite contar con estimaciones tempranas de variables que generalmente se conocen con cierto desface o retraso considerable.

El nowcasting se vale de información de variables relacionadas con la variable objetivo y de las cuales si contamos con información de manera oportuna.

Esta tecnica tiene aplicaciones en diversos ámbitos como:

  • Economía: Calculo PIB, inflación, desempleo.

  • Fianzas: Precios de acciones, bonos, divisas.

  • Salud: Propragación de enfermedades.

  • Medio ambiente: Clima, calidad del aire,

  • Retail: Demanda del consumidor.

Choi y Varian (2009, 2012) demostraron que los datos de búsqueda de Google podrían usarse como un señal externa para un modelo de predicción inmediata, pero sus métodos requerien seleccionar cuidadosamente el conjunto de predictores a utilizar.


Objetivo

El objetivo del presentre trabajo es analizar el comportamiento de la serie temporal de datos con cifras de ventas de casas nuevas y ajustar un modelo bayesiano de series de tiempo (BSTS) para hacer predicciones inmediatas y selección de variables explicativas, el cual contemple los siguientes aspectos:

  • Un módulo de series de tiempo, que capture la tendencia general y los patrones estacionales en los datos.

  • Un componente de regresión que permita la incorporación de datos exógenos (e.g. datos de Google Trends).

Metodología

Obtención y pre-procesamiento de datos

Los datos empleados para el análisis se obtuvieron de las estadísticas ,publicadas mensualmente por La Oficina del Censo de EE. UU. (The US Census Bureau) y el Departamento de Vivienda y Desarrollo Urbano de EE. UU. (The US Department of Housing and Urban Development), sobre el mercado de la vivienda al cierre de cada mes. Los datos incluyen cifras sobre ‘Casa nueva vendida y por Venta’ para el periodo que comprende desde enero del año 2004 hasta septiembre del 2009, se presentan de forma mensual y en unidades de miles.

Una vez obtenida la base anterior con 105 observaciones se emplearon las herramientas de google para obtener una base con términos de búsqueda de Google Trends (obtenidos de http://trends.google.com/correlate). Estos muestran la popularidad relativa de cada término de búsqueda entre todos los términos de búsqueda ingresados en Google. Lo que se hizo para lograr esto es usar la serie de tiempo de ventas de casas como input para Google Correlate y el output es otra serie con las 100 búsquedas de Google más correlacionadas. Finalmente se eliminan 31 variables del output que no hacen sentido, por ejemplo: “tahitian noni house”, “exhaust sound”; por lo que el modelo inicial tendrá sólo 69 regresores.

Las cifras en esta base obtenida con la herramienta de Google están estandarizadas por lo que se estandariza la nuestra variable objetivo y se unen las dos bases de acuerdo a la fecha . El resultados es una sola base con 105 renglones (observaciones) y 71 columnas (fecha, ventas de casas y 69 covariables).

data_housing <- read.csv("./data/HSN1FNSA_correlate.csv")
head(data_housing |> select(1:10))
##         Date HSN1FNSA X80.20.mortgage appreciation.rate home.appreciation
## 1 2004/01/01    0.942           1.080             0.380             1.198
## 2 2004/02/01    1.336           0.930             1.139             1.096
## 3 2004/03/01    1.973           2.091             1.734             1.520
## 4 2004/04/01    1.548           1.534             1.072             1.419
## 5 2004/05/01    1.730           1.595             1.524             1.584
## 6 2004/06/01    1.427           1.350             1.518             1.467
##   help.u.sell planned.community new.home.builder appreciation.rates city.land
## 1       0.724             1.682            1.371              1.173     0.987
## 2       0.813             1.770            1.460              1.119     1.472
## 3       1.523             1.653            1.894              1.356     1.535
## 4       1.603             1.341            1.619              1.152     1.509
## 5       1.502             1.877            1.648              1.181     1.902
## 6       1.385             1.603            1.529              1.304     1.748


## Serie de Tiempo

Visualización de viviendas nuevas vendidas

Tendencia

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Gráficando las ventas de casas nuevas, observamos que alcanzaron su punto máximo en marzo de 2005.

  • Su principal característica es la tendencia decreciente.

  • Esta serie también tiene algo de estacionalidad mensual.

Estacionalidad (mensual)

Podemos observar un patrón fuerte de estacionalidad en los primeros tres años. Con el mes de Marzo con un pico de ventas importante, una disminución el resto del año con un ligero repunte en el mes de Octubre. A partir de Abril 2007 se observa un comportamiento atípico relacionado a la crisis sufrida en el mercado americano. Para los años 2011 y 2012 podemos observar nuevamente una ligera estacionalidad mensual.

Autocorrelación

Esta ACF muestra correlaciones totales, capturando tanto el efecto de la tendencia, como la estacionalidad. Al tener una tendencia tan fuerte como veíamos en las gráficas anteriores, hace sentido que la serie tenga muy alta correlación con los meses inmediatos anteriores, en otras palabras el aspecto de la tendencia esta dominando nuestra gráfica ACF.

Podríamos también preguntarnos por correlaciones parciales, es decir podriamos interesarnos por ver el efecto que tienen la ventas del mes de junio de 2012 en las ventas de septiembre de 2012 cotrolando por las ventas de julio y agosto de 2012.

Es decir, buscamos la correlación directa, bloqueando la correlación que ocurre a través de otros meses intermedios.

Cuando condicionamos meses intermedios en nuestra serie, prácticamente todas las correlaciones parciales se encuentran en los intervalos del 95% (estos intervalos consideran la hipotesis de que los valores de la serie no tienen autocorrelación). Estos datos podrían corresponder a un modelo subyacente autoregresivo de orden 1.

Modelo bayesiano estructural de series de tiempo

Los modelos bayesianos estructurales de series de tiempo son modelos utilizados para la selección de covariables, forecasting, nowcasting, la inferencia del impacto causal y otras aplicaciones. La representación estado-espacio considerada para el modelo que se busca implementar es:

Donde \(\mu_t\) representa la tendencia, \(\tau_t\) la estacionalidad y \(\beta^Tx_t\) el componente de regresión. El componente de la tendencia es parecida a la ecuación de nivel pero con un término adicional \(\delta_t\) que denota la cantidad extra de \(\mu\) cuando damos un paso \(t \rightarrow t+1\) y puede interpretarse como la pendiente de la tendencia lineal local multiplicada por \(\Delta t\) que siempre es igual a 1. En general un modelo con tendencia lineal local es un modelo mejor que el modelo de nivel local si se cree que la serie de tiempo tiene una tendencia en una dirección particular y desea que los pronósticos futuros reflejen un aumento (o disminución) continuo visto en observaciones recientes. Mientras que el modelo de nivel local basa los pronósticos en torno al valor promedio de las observaciones recientes, el modelo de tendencia lineal local también agrega pendientes ascendentes o descendentes recientes. POr este motivo, para este trabajo buscamos comparar estos dos modelos.

La mejor forma de comprender el componente estacional \(\tau_t\) es verlo como una regresión con variables estacionales dummies. En nuestro caso tenemos periodos mensuales, de manera que S=12. El modelo de estados estacional incluye las 12 variables dummies pero restringe sus coeficientes para que sumen cero. Para el modelo de desempleo, Scott y Varian describieron el ciclo anual con los datos semanales de reclamaciones iniciales utilizando un componente estacional con S=52. Sabemos que las semanas dividen de manera perfecta a los años, pero dado el pequeño número de años para los cuales hay datos de Google disponibles, se consideró que la discontinuidad estacional ocasional de un período es irrelevante.

Para lo antes mencionado, los modelos con enfoque bayesiano que buscamos implementar y comparar ,con la paquetería de BSTS en R, constan de tres técnicas principales:

  • Filtro Kalman. La técnica de descomposición de series temporales. En este paso, un investigador puede agregar diferentes variables de estado: tendencia, estacionalidad y otras.

  • Método de Spike-and-Slab. En este paso, se seleccionan los predictores de regresión más importantes.

  • Media del modelo bayesiano. Combinando los resultados y el cálculo de predicción.

Selección de Modelo

Estamos interesados en hacer varios modelos y elejir el que mejor ajuste, uno de estos modelos debe incluir las variables regresoras por lo que eliminamos los últimos 3 meses de nuestros datos que después queremos predecir con este modelo.

Se realizó un modelo inicial en donde solo se contempla la tendencia.

#para los estados con tendencia
new_data <- tail(data_housing, 3)
data_housing <- tail(data_housing, 102)
ss1 <- AddLocalLinearTrend(list(), data_housing$HSN1FNSA)
model1 <- bsts(data_housing$HSN1FNSA,
               state.specification = ss1,
               niter = 2000)
## =-=-=-=-= Iteration 0 Wed May 17 00:27:51 2023 =-=-=-=-=
## =-=-=-=-= Iteration 200 Wed May 17 00:27:51 2023 =-=-=-=-=
## =-=-=-=-= Iteration 400 Wed May 17 00:27:52 2023 =-=-=-=-=
## =-=-=-=-= Iteration 600 Wed May 17 00:27:52 2023 =-=-=-=-=
## =-=-=-=-= Iteration 800 Wed May 17 00:27:52 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1000 Wed May 17 00:27:53 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1200 Wed May 17 00:27:53 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1400 Wed May 17 00:27:53 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1600 Wed May 17 00:27:54 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1800 Wed May 17 00:27:54 2023 =-=-=-=-=
plot(model1)

plot(model1, "components")

pred1 <- predict(model1, horizon = 3)
plot(pred1, plot.original = 102)


En un segundo modelo nos interesa tomar en cuenta la estacionalidad con periodo fijo de 12 y es muy claro como hay mejoría de predicción, el intervalo de confianza se hace más pequeño para los últimos 3 meses.

ss2 <- AddSeasonal(ss1, data_housing$HSN1FNSA, nseasons = 12)
model2 <- bsts(data_housing$HSN1FNSA,
               state.specification = ss2,
               niter = 2000)
## =-=-=-=-= Iteration 0 Wed May 17 00:27:55 2023 =-=-=-=-=
## =-=-=-=-= Iteration 200 Wed May 17 00:27:56 2023 =-=-=-=-=
## =-=-=-=-= Iteration 400 Wed May 17 00:27:57 2023 =-=-=-=-=
## =-=-=-=-= Iteration 600 Wed May 17 00:27:57 2023 =-=-=-=-=
## =-=-=-=-= Iteration 800 Wed May 17 00:27:58 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1000 Wed May 17 00:27:58 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1200 Wed May 17 00:27:59 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1400 Wed May 17 00:28:00 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1600 Wed May 17 00:28:00 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1800 Wed May 17 00:28:01 2023 =-=-=-=-=
plot(model2)

plot(model2, "components")

pred2 <- predict(model2, horizon = 3)
plot(pred2, plot.original = 102)

Un tercer modelo contempla un componente de autocorelación 1 basado en las gráficas de autocorrelación y autocorrelación parcial previamente mostradas. Vemos en las siguientes gráficas que el intervalo de estimación para el componente trend se hace muy grande, lo cual de incio no es un buen indicador.

ss3 <- AddAr(ss2, data_housing$HSN1FNSA, lags = 1)
model3 <- bsts(data_housing$HSN1FNSA,
               state.specification = ss3,
               niter = 2000, seed=12)
## =-=-=-=-= Iteration 0 Wed May 17 00:28:03 2023 =-=-=-=-=
## =-=-=-=-= Iteration 200 Wed May 17 00:28:03 2023 =-=-=-=-=
## =-=-=-=-= Iteration 400 Wed May 17 00:28:04 2023 =-=-=-=-=
## =-=-=-=-= Iteration 600 Wed May 17 00:28:05 2023 =-=-=-=-=
## =-=-=-=-= Iteration 800 Wed May 17 00:28:06 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1000 Wed May 17 00:28:06 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1200 Wed May 17 00:28:07 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1400 Wed May 17 00:28:08 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1600 Wed May 17 00:28:09 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1800 Wed May 17 00:28:09 2023 =-=-=-=-=
plot(model3)

plot(model3, "components")

pred3 <- predict(model3, horizon = 3)
plot(pred3, plot.original = 102)

Comparamos los tres modelos anteriores con la siguiente gráfica y vemos que los errores acumulados para el modelo que incluye la componente de autoregresón es ligeramente mayor a la de los dos primeros modelos. Esto nos ayuda a tomar la desición de descartar el tercer modelo. Pareciera tambien que el modelo que incluye la estacionalidad tiene errores acumulados ligramente mayor que el que solo toma la tendencia pero debido a que en la gráfica de este modelo se aprecia menor varianza en la predicción no lo desechamos.

CompareBstsModels(list("Modelo - trend" = model1,
                       "Modelo - trend + season" = model2,
                       "Modelo - trend + season + AR1" = model3),
                  colors = c("black", "red", "blue"))


Para definir si las variables obtenidas con las búsquedas de Google son útiles en la predicción inmediata que buscamos se realizó un un último modelo que considera el componente de regresión con 69 coeficientes para las covariables que queremos usar como regresoras. Este modelo solo toma en cuenta la estacionalidad y tendencia. El modelo también toma tamaño esperado inicial de 5, lo que corresponde al ajuste la probabilidad de inclusión a priori de 5/69 = 0.072.

#por default busca sólo tomar una variable
model4 <- bsts( HSN1FNSA~ .,
               state.specification = ss2,
               niter = 2000,
               data = select(data_housing, -1),
               expected.model.size = 5)
## =-=-=-=-= Iteration 0 Wed May 17 00:28:12 2023 =-=-=-=-=
## =-=-=-=-= Iteration 200 Wed May 17 00:28:13 2023 =-=-=-=-=
## =-=-=-=-= Iteration 400 Wed May 17 00:28:13 2023 =-=-=-=-=
## =-=-=-=-= Iteration 600 Wed May 17 00:28:14 2023 =-=-=-=-=
## =-=-=-=-= Iteration 800 Wed May 17 00:28:15 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1000 Wed May 17 00:28:15 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1200 Wed May 17 00:28:16 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1400 Wed May 17 00:28:17 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1600 Wed May 17 00:28:17 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1800 Wed May 17 00:28:18 2023 =-=-=-=-=
plot(model4)

plot(model4, "components")

# model3
pred4 <- predict(model4, horizon = 3, newdata = new_data)
plot(pred4, plot.original = 102)


Nos interesa ver la selección de variables

plot(model4, "coefficients", inc=.15, main = "Top de variables seleccionadas para el modelo 4")


Con la gráfica anterior podemos pensar en un modelo que de tamaño esperado por lo que nos interesa crear ahora un último modelo con tamaño esperado igual a dos:

#por default busca sólo tomar una variable
model5 <- bsts( HSN1FNSA~ .,
               state.specification = ss1,
               niter = 2000,
               data = select(data_housing, -1),
               expected.model.size = 2)
## =-=-=-=-= Iteration 0 Wed May 17 00:28:21 2023 =-=-=-=-=
## =-=-=-=-= Iteration 200 Wed May 17 00:28:21 2023 =-=-=-=-=
## =-=-=-=-= Iteration 400 Wed May 17 00:28:22 2023 =-=-=-=-=
## =-=-=-=-= Iteration 600 Wed May 17 00:28:22 2023 =-=-=-=-=
## =-=-=-=-= Iteration 800 Wed May 17 00:28:23 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1000 Wed May 17 00:28:23 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1200 Wed May 17 00:28:23 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1400 Wed May 17 00:28:24 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1600 Wed May 17 00:28:24 2023 =-=-=-=-=
## =-=-=-=-= Iteration 1800 Wed May 17 00:28:25 2023 =-=-=-=-=
plot(model5)

plot(model5, "components")

# model3
pred5 <- predict(model5, horizon = 3, newdata = new_data)
plot(pred5, plot.original = 102)

Para saber si realmente hay una mejora usando los datos de Google comparamos los tres modelos, que no descratamos, con la gráfica de errores acumulados de predicción a un paso y vemos que el modelo de tamaño dos que incluye las variables explicativas tiene mejor desempeño.

CompareBstsModels(list("Trend" = model1,
                       "Trend + Season" = model2,
                       "Trend + Season + 2Vars" = model5),
                  colors = c("black", "red", "blue"))


En esta gráfica se puede ver que para el año 2005 que es donde empieza una tendecia a la baja más notable, y que corresponde al periodo 28, la diferencia de los errores acumulados se va haciendo más grande entre los modelos que no incluyen variables contra los que si.

Si quisieramos verificar si los residuos para este modelo se pueden ver como ruido blanco podemos verlo con la siguiente gráfica.

pred_errors_tbl <- bsts.prediction.errors(model5)$in.sample |> 
  t() |> as_tibble() |>
  mutate(t = 1:102) |> 
  pivot_longer(-c(t), names_to = "sim", values_to = "valor") |>
  group_by(t) |> 
  summarise(valor = mean(valor)) |> 
  as_tsibble(index = t)

ACF(pred_errors_tbl, valor, lag_max = 20) |> 
  autoplot() + ylim(c(-1,1))



Conclusión

El halcón murcielaguero presenta una amplia distribución en México, la cual está altamente correlacionada con la idoneidad ambiental para la especie. Sus poblaciones más abundantes y extensas se encuentran en el sureste del país, ## Referencias

  • Banbura, Marta and Giannone, Domenico and Reichlin, Lucrezia, Nowcasting (Noviembre 30, 2010). ECB Working Paper No. 1275, Disponible en SSRN: https://ssrn.com/abstract=1717887 or http://dx.doi.org/10.2139/ssrn.1717887
  • Choi, H. and Varian, H. (2009). Predicting the present with Google Trends. Tech. rep., Google.
  • Choi, H. and Varian, H. (2012). Predicting the present with Google Trends. Economic Record 88, 2–9.
  • U.S. Census Bureau and U.S. Department of Housing and Urban Development, New One Family Houses Sold: United States [HSN1FNSA], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/HSN1FNSA